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ABSTRACT 


We  have  developed  a  3-D  velocity  model  for  the  crust  and  upper  mantle  in  the  India-Pakistan  region 
(WINPAK3D)  to  improve  regional  event  location.  Results  of  extensive  testing  demonstrate  that  the  3-D 
model  improves  location  accuracy  for  this  region,  specifically  for  the  case  of  small,  regionally  recorded 
events,  for  which  teleseismic  data  may  not  be  available.  The  initial  model  for  the  region  was  developed  by 
integrating  the  results  of  more  than  sixty  previous  studies  related  to  crustal  and  upper  mantle  velocity 
structure.  During  the  final  year  of  this  effort,  we  have  developed  a  3-D  joint  velocity  tomography/event 
relocation  to  improve  this  preliminary  velocity  model.  The  algorithm  applies  an  iterative,  conjugate- 
gradients  technique  to  minimize  the  misfit  between  calculated  and  observed  travel  times  from  multiple 
stations  and  events,  subject  to  smoothness  and  geologic  constraints.  Events  are  relocated  using  the  grid 
search  location  method  of  Rodi  and  Toksdz  (2000)  after  each  update  of  the  3-D  velocity  model.  Travel 
times  and  their  sensitivities  to  the  velocity  structure  are  computed  with  an  extension  of  the  3-D  Podvin- 
Lecomte  (1991)  method.  We  are  applying  our  tomography  algorithm  to  a  suite  of  580  events  containing 
over  7,800  arrivals  from  the  Engdahl  et  al.  (1998)  database. 

A  critical  aspect  of  model  development  is  validation,  which  we  have  addressed  in  three  ways:  (1)  cross- 
validation  analysis  for  a  variety  of  events,  (2)  comparison  of  model-determined  hypocenters  with  known 
event  location,  and  (3)  comparison  of  model-derived  and  empirically  derived  source-specific  station 
corrections  (SSSC's)  generated  for  the  International  Monitoring  System  (IMS)  auxiliary  seismic  station 
located  at  Nilore,  Pakistan.  The  3-D  model  provides  improvement  in  regional  location  compared  to  a 
global  1-D  model.  For  example,  an  event  with  a  well-known  location  had  an  epicenter  mislocation  of  only 
6. 1  km  for  the  solution  determined  using  the  3-D  model,  compared  with  a  mislocation  of  15.4  km  for  the 
IASPEI91  model  (Kennett  and  Engdahl,  1991).  These  and  other  results  demonstrate  that  3-D  models  are  a 
prerequisite  for  achieving  improved  location  accuracies  in  areas  with  complex  crustal  and  upper  mantle 
structure. 

KEY  WORDS:  3-D  velocity  model,  India,  Pakistan,  tomography,  improved  event  location 


OBJECTIVE 


Improving  seismic  event  location  is  particularly  important  for  monitoring  the  Comprehensive  Nuclear-Test- 
Ban  Treaty  (CTBT).  The  treaty  requires  a  location  accuracy  of  better  than  1000  km^  due  to  the  restrictions 
defined  both  for  on-site-inspections  (OSI)  and  for  using  location  as  a  discrimination  tool.  For  small  events 
(mt,  <  4.0),  meeting  this  requirement  is  problematic  because  of  limited  recordings  and  complicated  crustal 
structure  at  regional  distances  not  accounted  for  in  the  standard  global  1-D  models  such  as  the  IASPEI91 
model  (Kennett  and  Engdahl,  1991).  Precise  regional  location  of  seismic  events  in  the  context  of  the  CTBT 
requires  a  velocity  model  that  accurately  represents  the  real  Earth  structure,  as  systematic  biases  caused  by 
unmodeled  Earth  structure  are  known  to  play  an  important  role  in  earthquake  location  errors  (e.g.,  Douglas, 
1967;  Dewey,  1972;  Engdahl  and  Lee,  1976;  Jordan  and  Sverdrup,  1981;  Pavlis,  1992).  Regional  3-D 
models  that  better  represent  the  true  Earth  structure  can  be  used  to  compute  accurate  travel  times  of 
regional  seismic  phases  such  as  Pn,  Pg,  Sn,  and  Lg.  Travel  times  calculated  using  3-D  models  can  then  be 
used  to  develop  source  specific  station  corrections  (SSSC’s)  that  can  be  implemented  by  monitoring 
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organizations  to  provide  improved  regional  event  loeations.  Our  goal  is  to  develop  a  regional  model  of  the 
crust  and  upper  mantle  for  the  India-Pakistan  region  that  will  improve  event  location  accuracy  and  provide 
a  set  of  improved  hypocenter  estimates  based  on  that  model.  We  are  building  an  accurate  3-D  model  for 
the  India-Pakistan  region  by  updating  a  preliminary  model  through  a  joint  velocity  tomography  and 
hypocenter  relocation  technique.  Here  we  describe  this  approach  in  detail  and  our  progress  towards 
developing  an  accurate  model  for  the  India-Pakistan  region. 

RESEARCH  ACCOMPLISHED 


Preliminary  Model 

We  have  developed  a  regional  3-D  velocity  model  of  the  crust  and  upper  mantle  structure,  the  Weston 
INdia/PAKistan  3-D  model  (WINPAK3D),  to  improve  seismic  event  locations  for  the  India-Pakistan 
region  (Figure  la).  The  region  of  interest  has  a  complex  tectonic  history  exhibited  in  the  diverse 
geometries  of  crustal  structures  across  the  area.  Using  the  available  geophysical  literature,  we  have 
compiled  a  detailed  3-D  velocity  model  for  the  India-Pakistan  region  (Figure  Ib).  We  synthesized  data 
from  approximately  sixty  published  references  relevant  to  the  velocity  structure,  geology,  and  tectonics 
throughout  the  region.  These  references  include  data  such  as  seismic  reflection  and  refraction  surveys  (i.e. 
DSS  profdes,  Pn  tomography,  Pnl  waveform  inversion),  interpretations  of  gravity  data,  surface  wave 
studies,  and  receiver  function  analyses.  Because  these  data  sources  vary  in  spatial  coverage,  resolution,  and 
the  number  of  constraints,  the  model  varies  in  a  similar  manner.  The  velocity  model  is  defined  on  a  grid  of 
one-degree  by  one-degree  blocks  and  5  km  depth  intervals  from  0  to  75  km.  We  have  appended  the 
IASPEI91  model  (Kennett  and  Engdahl,  1991)  to  the  base  of  the  preliminary  velocity  model,  beginning  at 
80  km  depth  and  extending  to  700  km  depth,  to  accommodate  raypaths  that  travel  into  the  upper  mantle. 
See  Johnson  and  Vincent  (2001)  for  additional  details  on  the  development  of  the  WINPAK3D  model. 


Figure  1.  (a)  Regional  setting  of  the  WINPAK3D  velocity  model  showing  the  location  of  seismic  station 
NIL  at  Nilore,  Pakistan  (future  site  of  IMS  station  PRPK)  as  the  black  triangle  and  locations  of  the 
India  and  Pakistan  nuclear  test  sites  as  black  stars,  (b)  Slices  of  the  velocity  model  at  selected 
depths  show  the  variability  of  velocity  model  resolution,  as  well  as  the  large  range  of  crustal 
thickness  across  this  complex  tectonic  region. 
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Evaluating  the  Initial  Model 

The  initial  WINPAK3D  velocity  model  was  evaluated  to  ensure  a  reasonable  model  as  a  starting  point  for 
tomographic  inversion.  Testing  and  validation  of  a  velocity  model  is  a  difficult  task  for  several  reasons: 

(1)  the  paucity  of  ground  truth  data,  Le.,  explosions  and  earthquakes  with  known  or  well  constrained 
location,  as  defined  by  Yang  et  al.  (2000),  (2)  the  pitfalls  of  using  root  mean  square  (rms)  residuals  to 
indicate  improved  location  accuracy  ,  and  (3)  the  need  to  avoid  circularity  issues  by  ensuring  that  data  used 
to  generate  the  model  are  omitted  from  tests  to  validate  the  model.  We  evaluate  our  regional  model's 
validity  by  comparing  the  hypocenter  from  a  well-located  event,  recorded  at  a  number  of  regional  seismic 
stations,  with  model  derived  hypocenters.  In  addition,  we  implement  a  robust  statistical  procedure  to 
further  examine  the  validity  of  our  regional  model.  We  compare  the  performance  of  the  3-D  model  to  the 
IASPEI91  model  to  determine  the  improvement  over  the  global  1-D  model  for  this  region.  Since  the  model 
was  under  development  at  the  time  of  testing,  the  southern  one-third  of  the  current  model  was  not  included 
in  the  following  model  evaluation  tests. 

We  begin  by  comparing  our  regional  3-D  velocity  model  with  the  IASPEI91  model  using  the  cross 
validation  technique,  a  statistical  method  of  evaluating  a  model.  Cross  validation  obtains  nearly  unbiased 
estimators  of  prediction  error  in  complicated  problems  (Efron  and  Gong,  1983).  We  performed  this 
analysis  for  6  events  (Figure  2),  chosen  because  they  are  well  recorded  and  their  epicenters  are  distributed 
across  the  model  region.  Applying  the  cross  validation  approach  to  the  location  problem,  we  repeatedly 
solved  for  the  hypocenter  of  each  event  using  unique  subsets  of  5  stations  (simulating  a  sparse  array) 
randomly  selected  from  the  set  of  regional  arrivals.  For  each  realization,  we  predicted  the  travel  times  to 
the  omitted  stations  based  on  the  hypocenter  determination  derived  from  the  5  stations.  We  calculated  the 
average  residual  of  predicted  travel  times  at  each  station  for  all  realizations  and  then  evaluated  the  model  fit 
in  terms  of  the  rms  error  for  all  stations.  If  the  devised  model  is  a  physically  relevant  one,  then  it  will  yield 
accurate  predictions  for  the  omitted  data. 


•  Event  1 

•  Event  2 

•  Event  3 

•  Event  4 

•  Event  5 

•  Event  6 


Figure  2.  Events  used  for  the  cross  validation  analysis  to  test  travel  time  prediction  capabilities  of  the 

WINPAK3D  model  as  compared  with  IASPEI91.  Raypaths  to  the  regional  stations  that  recorded 
the  events  are  shown. 

We  performed  this  analysis  for  both  the  WINPAK3D  and  IASPEI91  models  and  calculated  the  ms  error 
for  the  travel  time  predictions  made  from  each  model  (see  Table  1).  WINPAK3D  achieves  a  smaller 
prediction  error  for  the  India/Pakistan  region  for  5  of  the  6  events,  with  a  large  range  of  improvement  in 
travel  time  prediction  accuracy  (22-83%)  over  IASPEI91.  For  event  one,  the  prediction  was  the  same  for 
both  WINPAK3D  and  IASPEI91.  Although  the  ms  prediction  results  demonstrate  that  WINPAK3D 
improves  travel  time  prediction  accuracy  over  IASPEI91  for  events  and  raypaths  throughout  the  model 
region,  the  amount  of  improvement  that  the  3-D  model  provides  over  IASPEI91  is  quite  variable.  This  is 
likely  due  to  the  varying  resolution  of  the  model  discussed  in  the  previous  section.  These  results  suggest 
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that  WINPAK3D  will  serve  as  a  good  starting  model  for  the  tomography  inversion,  yet  eertain  regions  of 
the  model  will  be  subjeet  to  bigger  updates  than  others. 

Table  1.  Cross  Validation  Results:  RMS  Predietion  Error 

Model  Prediction  Error  (seconds) 


Event 

Stations 

Realizations 

WINPAK3D 

IASPEI91 

Improvement 

1 

5/14 

2002 

2.1 

2.1 

0% 

2 

5/13 

1287 

6.4 

8.7 

26% 

3 

5/17 

5000 

6.8 

8.7 

22% 

4 

5/18 

5000 

5.7 

7.4 

23% 

5 

5/14 

2002 

3.9 

6.9 

43% 

6 

5/23 

5000 

1.3 

7.6 

83% 

Event  6  (mt,  5.2,  USGS)  in  the  eross  validation  analysis  oeeurred  on  14  February  1977  in  the  region  near 
Nilore,  Pakistan  (Figure  2).  This  event  was  well  reeorded  due  to  the  presenee  of  two  temporary  networks, 
Tarbela  and  Chasma,  and  a  subarray  of  five  stations  eentered  on  Nilore,  Pakistan  (NIL),  all  loeated  in  elose 
proximity  to  the  event .  Based  solely  on  data  from  these  three  networks,  Seeber  and  Armbmster  (1979) 
estimated  the  hypoeenter  of  the  event  and  eondueted  a  detailed  study  of  the  aftershoek  sequenee.  Using 
teleseismie  as  well  as  regional  data  (not  ineluding  the  temporary  networks),  both  the  International 
Seismologieal  Centre  (ISC)  and  the  U.S.  Geologieal  Survey  (USGS)  loeated  the  event  within 
approximately  5  km  of  the  epieenter  reported  by  Seeber  and  Armbmster  (S-A).  The  hypoeenters  of  the 
main  shoek  (depth  =  14.46  km;  John  Armbmster,  personal  eommunieation,  2000)  and  the  50  aeeurately 
loeated  aftershoeks  determined  by  S-A  indieate  a  mpture  surfaee  between  12  and  18  km  depth.  Given  the 
high  loeation  aeeuraey  estimated  for  this  event  as  well  as  the  large  number  of  regional  arrivals,  we  ehose 
this  event  for  a  test  of  our  model's  loeation  eapabilities. 

The  S-A  loeation  (1979)  for  this  event  is  one  of  the  better  eonstrained  loeations  in  the  region  due  both  to 
the  loeal  network  eoverage  and  to  the  detailed  aftershoek  study.  Though  there  have  been  several  nuelear 
tests  in  the  region,  elassified  as  ground  tmth  (Yang  et  ah,  2000)  with  known  aeeuraeies  of  0  km  and  1  km 
(GTO  and  GTl),  the  regional  data  availability  within  our  model  area  is  limited.  Arrival  times  for  the  14 
Febmary  1977  event  were  reported  to  the  ISC  by  23  seismie  stations  within  our  model  region.  This  event 
is  additionally  important  due  to  its  proximity  to  the  seismie  station  at  Nilore,  Pakistan  (NIL),  the  future  site 
of  the  International  Monitoring  System  (IMS)  station  PRPK. 


Table  2.  Location  accuracy  results  of  the  14  Febmary  1977  Pakistan  event 


Reference 

RMS 

Residnal  (s) 

Eat. 

Lon. 

Depth  (km) 

Mislocation  (km) 
Epicenter  Depth 

Seeber  and  Armbmster 

0.05 

33.625 

73.208 

14.46 

- 

- 

WINPAK3D 

1.06 

33.57 

73.20 

11.1 

6.1 

-3.3 

IASPEI91 

1.47 

33.51 

73.11 

8.9 

15.4 

-5.5 

NEIC 

- 

33.595 

73.253 

33 

5.0 

18.5 

ISC 

1.25 

33.597 

73.267 

27 

5.2 

12.5 

We  believe  that  the  velocity  model  should  be  validated  using  many  source-station  pairs  prior  to  assessing 
model  performance  with  sparse  data.  Using  the  23  regional  arrivals,  we  calculated  the  hypoeenters 
predicted  by  both  the  WINPAK3D  and  IASPEI91  models  (Table  2)  for  this  event.  The  epicenter 
determined  by  WINPAK3D  is  6.1  km  from  the  "best"  event  location  as  determined  by  S-A  (1979),  while 
the  IASPEI91  epicenter  is  approximately  15  km  from  the  event.  The  depth  is  well  determined  for  both 
models,  due  to  extensive  near-regional  coverage  when  data  from  all  stations  are  used.  However,  the  cross 
validation  results  suggest  that,  in  stmcturally  complex  regions,  global  models  will  be  unacceptable  for 
accurate  location  of  small  events  not  recorded  teleseismically. 

Therefore,  we  compare  the  locations  for  the  event  as  determined  by  both  the  WINPAK3D  and  IASPEI91 
models  using  only  four  stations,  serving  as  IMS  station  surrogates.  The  four  surrogate  stations  and  their 
IMS  counterparts  are  shown  in  Figure  3.  The  location  determined  using  WINPAK3D  and  only  arrivals 
from  VAN,  FRU,  MNL,  and  TGI  (Figure  4a,  Table  3)  is  remarkably  accurate  (4.6  km).  The  IASPEI91 
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solution  is  within  14  km  of  the  S-A  epicenter.  However,  the  surface  projection  of  the  3-D  95%  confidence 
region  determined  by  the  Monte  Carlo  simulation  algorithm  of  Rodi  and  Toksoz  (2000),  encompasses  the 
S-A  epicenter  for  the  WINPAK3D  model,  but  not  for  the  IASPEI91  model  (Figure  5a).  Furthermore,  the 
epicenter  calculated  using  the  3-D  model  and  only  four  regional  data  is  as  accurate  as  the  epicenters 
produced  by  both  the  ISC  and  the  USGS  (Figure  4a),  which  included  teleseismic  as  well  as  regional  data. 


Figure  3.  Location  of  the  14  February  1977  event  (small  circle  plotted  on  top  of  PRPK)  and  the  locations 
of  regional  stations  VAN,  FRU,  MNL,  and  TGI  (black  triangles)  used  as  surrogates  for  IMS 
stations  (red,  inverted  triangles)  for  the  sparse  array  locations.  Blue  triangles  (LAH,  BHK,  DDI, 
and  NDI)  show  stations  used  in  the  sparse  array  depth  accuracy  study. 

We  used  the  same  four  stations  to  determine  the  event  location  while  fixing  the  depth  to  15  km,  the 
approximate  depth  of  the  S-A  (1979)  solution.  With  the  hypocenter  fixed  at  the  approximate  event  depth, 
epicenter  mislocation  decreases  to  2.9  km  for  the  WINPAK3D  model  while  it  increases  to  17  km  for  the 
IASPEI91  model.  The  confidence  regions  for  fixed  depth  solutions  are  smaller,  and  again,  the 
WINPAK3D  confidence  region  encompasses  the  S-A  location  while  the  IASPEI91  confidence  region  does 
not.  Interestingly,  the  solutions  determined  by  both  models,  for  both  the  fixed  depth  and  free  depth  cases, 
lie  within  the  1000  km^  region  (17.8  km  radius)  required  for  CTBT  monitoring  protocol.  We  believe  this  is 
due  to  the  close  proximity  of  MNL  to  the  event,  helping  to  constrain  the  event  depth. 

Next,  we  next  compare  depth  accuracy  between  WINPAK3D  and  IASPEI91  for  the  4  station  problem 
where  different  stations  were  substituted  for  MNL,  the  station  closest  to  the  epicenter  (Figure  4).  Given 
that  the  event  is  known  to  be  a  shallow  crustal  event,  the  only  constraint  imposed  was  to  limit  depth  to  less 
than  100  km  depth.  The  results  for  the  4  station  problem  where  depth  is  not  controlled  by  a  local  station 
are  shown  in  Table  3.  The  3-D  model  provides  much  better  constraint  on  event  depth  than  the  IASPEI91 
model.  WINPAK3D  determined  the  depth  of  the  event  within  10  km  of  the  S-A  (1979)  depth  for  all  4 
variations  of  station  coverage.  However,  all  solutions  calculated  with  the  IASPEI91  model  placed  the 
event  at  the  maximum  allowed  depth  of  100  km.  Epicenter  mislocation  was  smaller  for  the  WINPAK3D 
model  as  well,  and  in  one  case  epicenter  mislocation  for  the  IASPEI91  solution  exceeded  the  17.84  km 
radius  required  for  the  1000  km^  OSI  region.  More  importantly,  the  IASPEI91  based  solutions  could 
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falsely  rule  out  this  shallow  event  based  on  the  erroneously  deep  hypoeenters.  Therefore,  for  this  event,  the 
WINPAK3D  model  has  been  shown  to  be  extremely  important  for  aehieving  the  loeation  aeeuraey  required 
in  nuelear  monitoring. 


Table  3.  Loeation  Aeeuraey  for  Sparse  Arrays 


WINPAK3D  Mislocation  (km) 

IASPEI91  Mislocation  (km) 

Stations 

Epicenter 

Depth 

Epicenter 

Depth 

IMS  Surrogates 

4.6 

-14.5 

13.5 

-14.5 

TGI,  VAN,  FRU,  LAH 

8.7 

0.5 

20.5 

85.5 

TGI,  VAN,  FRU,  BHK 

7.1 

-9.4 

10.5 

85.5 

TGI,  VAN,  FRU,  DDI 

1.9 

-9.4 

10.5 

85.5 

TGI,  VAN,  FRU,  NDI 

0.8 

0.5 

11.9 

85.5 

The  International  Data  Center  (IDC)  eurrently  uses  the  IASPEI91  global  travel  time  eurves  as  the  default 
for  event  loeation.  When  the  global  model  is  insuffieient  for  eharaeterizing  regional  geology,  the  IDC 
sometimes  applies  path  eorreetions  (Yang  et  al,  1999)  in  the  form  of  SSSC’s  that  are  a  funetion  of  souree 
loeation  for  any  station  and  phase.  We  derived  SSSC’s  based  on  the  WINPAK3D  veloeity  model  for  the 
India-Pakistan  region  as  deseribed  in  Reiter  et  al.  (2001).  To  aeeount  for  the  event  depth,  we  ealeulated  the 
antieipated  SSSC  for  station  NIL  based  on  the  WINPAK3D  model  for  an  event  at  15  km  depth. 


Figure  4.  Loeation  of  14  February  1977  event  as  determined  by  the  WINPAK3D  and  IASPEI91  models 
using  only  arrivals  from  four  stations  (VAN,  TGI,  FRU,  MNL)  serving  as  IMS  surrogates  when 
(a)  solving  for  hypoeenter  and  (b)  solving  for  epieenter  only,  with  depth  fixed  to  15  km.  The 
epieenter  given  by  Seeber  and  Armbruster  (1979)  is  indieated  by  the  star.  ISC  and  NEIC  solutions 
determined  using  both  regional  and  teleseismie  data  are  shown  by  the  open  eireles. 

The  14  February  1977  event  provides  a  unique  reeiproeity  test  of  the  3-D  veloeity  model,  sinee  it  oeeurred 
very  near  the  Nilore,  Pakistan  seismie  station  (NIL/PRPK).  The  residuals  with  respeet  to  IASPEI91  at  eaeh 
station  serve  as  empirieal  SSSC’s.  The  differenees  between  the  empirieal  and  model-based  SSSC’s  at  eaeh 
station  are  less  than  1  seeond  for  nearly  half  of  the  23  stations  and  less  than  3  seeonds  for  the  remaining 
stations.  Sinee  the  SSSC  values  are  as  large  as  +!-  8  seeonds  in  some  parts  of  the  region,  diserepaneies 
between  the  SSSC's  of  around  1  seeond  are  relatively  small  and  show  that  the  model-based  and  empirieal 
SSSC’s  for  NIL/PRPK  are  in  very  good  agreement. 

Tomography 

A  subset  of  580  events  were  seleeted  from  the  Engdahl  et  al.  (1998)  database  of  well  loeated  earthquakes 
and  explosions.  The  seleetion  eriteria  and  proeessing  of  events  in  the  Engdahl  et  al.  (1998)  data  set  results 
in  loeation  aeeuraey  of  15  km  or  better  in  most  eontinental  areas,  aeeording  to  the  lASPEI  Working  Group 
on  Referenee  Events  (http://lemond.eolorado.edu/~eopgte/).  The  subset  of  events  were  ehosen  in  order  to 
provide  the  best  spatial  distribution  of  events  aeross  the  region,  both  in  latitude-longitude  and  in  depth. 
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Additionally,  only  events  with  6  or  more  regional  arrivals  were  selected  from  the  database  to  insure 
sufficient  data  for  hypocenter  relocation.  These  events  were  recorded  at  78  stations  within  our  model 
region,  amounting  to  over  7,800  arrivals  (Figure  5).  Travel  times  from  these  events  will  be  used  in  an 
iterative  velocity  model  tomography  and  hypocenter  relocation  procedure  (described  below)  to  update  the 
preliminary  model. 


Figure  5.  Ray  coverage  for  tomography  inversion  provided  by  the  580  events  (circles)  recorded  at  78 
stations  (triangles).  The  region  outlined  by  the  black  rectangle  is  expected  to  have  the  best 
resolution  in  the  tomographic  update.  Events  and  stations  beyond  this  zone  were  included  to 
provide  optimal  ray  coverage  in  the  primary  model  region.  Ray  coverage  with  depth  is  also  very 
good. 

We  have  re-parameterized  WINPAK3D  to  reduce  the  number  of  model  parameters  to  be  solved  for  in  the 
inversion.  Our  new  parameterization  is  given  in  terms  of  a  velocity  vs.  depth  profile  at  each  point  on  a 
geographic  grid  sampled  unifonuly  in  latitude  and  longitude.  At  each  geographic  grid-point,  the  velocity 
profile  is  given  as  velocity/depth  pairs  at  nodes  ranging  from  sea  level  to  a  depth  of  760  km. 

Discontinuities  in  velocity  are  allowed  at  the  ocean  bottom,  Moho  and  the  major  mantle  discontinuities 
(410  and  660  km  in  the  IASP91  model). 

Our  initial  velocity  model  will  be  updated  to  fit  the  arrival  times  data  from  our  database  of  580  events.  To 
accomplish  this,  we  are  developing  an  inversion  algorithm  that  performs  joint  velocity  tomography  and 
multiple-event  location.  The  algorithm  will  relocate  the  earthquakes  in  conjunction  with  revising  the 
parameters  of  the  velocity  model.  Initially,  the  model  update  will  be  restricted  to  the  Pn  velocity  where  we 
expect  the  greatest  changes.  Subsequent  inversions  will  consider  updates  to  the  crustal  thickness  and  other 
parameters. 

Our  inversion  approach  is  formulated  as  follows.  The  unknowns  are  a  vector  m  containing  the  velocity 
model  parameters  to  be  estimated  (e.g.  Pn  velocity  at  each  point  of  the  latitude-longitude  grid),  and  the 
hypocenters  and  origin  times  of  several  events:  (xy,  t,),  j  =  1,...,  M.  The  data  are  arrival  times,  dy,  from 
each  event  to  a  subset  of  stations  indexed  as  i  =  1,...,  N.  The  data  and  unknowns  are  related  by 

dy  =  tj+Ti{^j-,m)  +  eij 
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where  ey  is  the  error  in  dy  and  T,  is  a  function  that  computes  traveltime  to  station  i  from  an  event 
hypocenter  Xy.  This  function  depends  on  the  model  parameter  vector  m.  Our  joint  inversion  criterion  is  to 
minimize  an  objective  function  of  the  form 

=  -t.  +7;.(x^.;m)|"  /4  +  r|Lm|" 

with  respect  to  all  the  unknowns.  Here,  <5y  is  the  standard  deviation  of  ey.  The  second  term  of  v|/  imposes  a 
smoothness  constraint  on  the  velocity  model,  with  L  being  a  regularization  operator  and  x  a  regularization 
parameter.  The  operator  L  is  chosen  as  a  differencing  operator  with  the  effect  that  spatial  derivatives  of  the 
model  velocity  are  minimized.  The  parameter  r  determines  the  degree  of  model  smoothness. 

Our  algorithm  minimizes  vj/  using  a  combination  of  conjugate  gradients  and  grid-search  techniques.  The 
conjugate  gradients  method  is  used  to  update  m  iteratively  along  a  sequence  of  computed  search  directions. 
At  each  conjugate  gradients  step,  grid  search  is  used  to  minimize  v|/  with  respect  to  the  event  origin 
parameters  with  m  fixed,  thus  updating  the  Xy  and  tj.  The  grid  search  for  a  given  event  is  performed  within 
a  specified  epicentral  radius  and  depth  range  from  its  initial  location,  allowing  us  to  handle  events  of 
varying  ground-truth  levels  (e.g.,  GTO,  GTS,  GT15).  Our  grid  search  algorithm  is  described  by  Rodi  and 
Toksdz  (2000). 

The  forward  model  for  this  inverse  problem  is  embodied  in  the  traveltime  functions  7).  For  fixed  x,  we 
evaluate  7)  (x;  m)  by  interpolating  a  traveltime  table  stored  on  a  3-D  hypocenter  grid.  The  grid  is  created 
by  applying  the  Podvin-Lecomte  (P-L)  finite-difference  traveltime  algorithm  (Podvin  and  Lecomte,  1991) 
to  the  earth  model  defined  by  m,  using  the  location  of  the  /th  station  as  the  "source"  in  the  calculation.  The 
forward  P-L  time  calculation  traces  the  wavefront  propagation  based  on  the  Huygens’  principle  in  the  finite 
difference  approximation.  The  model  is  discretized  with  an  equally  spaced  grid  comprised  of  constant 
velocity  cells.  Multiple  arrivals  (transmitted,  diffracted,  and  head  waves)  are  calculated  at  each  grid  node 
and  the  first  arrival  time  is  chosen.  The  time  t  at  the  current  node  is  a  function  of  the  times  t„  at  some  (3  or 
fewer)  of  the  neighboring  nodes  and  the  slowness,  s,  in  the  cell  traversed  by  the  wavefront  to  reach  this 
node.  That  is,  {tn,  ^). 

We  have  extended  the  P-L  algorithm  to  compute  the  sensitivities  of  traveltimes  to  cell  velocities,  as 
required  for  the  inversion.  When  performing  the  forward  calculation,  we  save  the  node  pattern  ("stencil") 
that  shows  which  of  the  neighboring  nodes  were  used  to  calculate  the  minimum  time  at  the  current  node. 
The  stencils  are  used  to  back-trace  from  any  node  of  the  grid  to  the  source.  The  ray-tracing  is  done  by 
identifying  all  of  the  grid  nodes  and  the  cells  (slownesses)  that  contribute  to  the  calculation  of  the  time  at 
the  receiver.  As  the  wavefront  propagates  away  from  the  source,  more  nodes  (and  cells)  are  involved  in  the 
travel  time  calculation  at  each  step.  After  the  midpoint  of  the  raypath,  the  propagation  region  narrows  until 
it  reaches  the  one  node  at  the  source  location. 

The  sensitivity  of  the  travel  time  to  the  slowness,  dtids,  is  calculated  at  each  grid  node  of  the  "ray"  for  the 
last  cell  traversed  by  the  wavefront  to  reach  that  node.  The  weight  of  each  neighboring  node  in  the 
calculation  of  the  time  at  the  current  node  ( dtidtn)  determines  the  weight  of  the  subsequent  node-to-source 
subpath  in  the  total  travel  time  calculation  for  the  ray.  The  sensitivities  along  each  subpath  are  weighted  by 
this  term. 

The  ray-tracing  algorithm  and  sensitivity  calculation  were  tested  by  comparing  the  travel  time  computed  as 
the  sum  of  the  weighted  sensitivity-slowness  products  to  the  forward  P-L  calculated  times.  For  rays  with 
10^  nodes,  the  difference  between  travel  times  calculated  by  these  two  methods  is  on  the  order  of  10’^  s 
when  the  calculation  is  done  in  single  precision.  Error  of  this  magnitude  for  the  number  of  nodes  in  the  ray 
can  be  accounted  for  by  the  level  of  precision  in  the  calculation,  demonstrating  that  the  ray-tracing 
technique  is  accurately  tracing  the  minimum  time  P-L  raypath. 

We  have  developed  the  necessary  algorithms  for  mapping  our  universal  velocity  model  to  a  Cartesian  block 
model  needed  by  the  P-L  algorithm,  and  for  mapping  the  3-D  Cartesian  traveltime  grids  to  geographic 
grids.  The  sensitivities  are  mapped  from  Cartesian  blocks  to  the  universal  model  parameterization. 
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We  note  that  our  joint  inversion  algorithm  is  fully  nonlinear  with  respect  to  both  the  velocity  model  and 
event  locations  since  traveltime  modeling  and  event  relocation  are  perfomied  for  each  update  of  m. 
However,  this  comes  at  a  high  computational  price  and  we  are  also  allowing  for  multiple  steps  of  the 
conjugate  gradients  iteration  before  the  P-L  modeling  is  repeated. 
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Figure  6.  Ray  sensitivities  for  station  DSH.  (a)  Location  of  station  DSH  (triangle)  and  raypaths  to  events 
(small  circles)  within  our  model  region  that  were  recorded  at  DSH.  Sensitivities  of  rays  to  model 
slices  at  (b)  50  km ,  (c)  60  km,  (d)  and  70  km  depth. 


Preliminary  Analysis  of  Ray  Sensitivities 

We  have  performed  a  test  of  our  ray  sensitivity  algorithm  by  computing  the  sensitivities  to  cell  velocities  of 
all  the  arrivals  from  our  event  database  that  were  recorded  at  station  DSH  (Figure  6a).  Selected  depth 
slices  of  the  sensitivities  for  DSH  are  shown  in  Figure  6b-6d.  The  rays  have  strong  sensitivity  to  the  model 
near  the  Moho  boundary  (presumably  in  the  Pn  velocity  depth  range).  This  suggests  that  updates  to  the  Pn 
velocities  as  well  as  the  Moho  depth  will  provide  the  most  significant  improvement  to  the  velocity  model. 

CONCLUSIONS  AND  RECOMMENDATIONS 


Validation  of  the  final,  tomographically  inverted,  regional  3-D  velocity  model  (WINPAK3D)  will  include 
more  extensive  application  of  the  methods  presented  within  this  paper  for  testing  of  the  preliminary  model. 
We  believe  that  all  three  techniques  are  important  to  the  validation  of  any  velocity  model.  Cross  validation 
provides  a  statistical  evaluation  of  the  model's  travel  time  prediction  capabilities.  Evaluation  of  location 
accuracy  for  a  number  of  ground  truth  events  is  a  very  important  aspect  of  model  validation,  as  it  will 
provide  estimates  of  the  model's  location  capabilities  for  location-based  discrimination.  Finally,  model 
versus  empirical  based  comparisons  such  as  the  comparison  of  the  WINPAK3D  SSSC  for  NIL  with  the 
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empirical  SSSC  derived  from  data  for  the  14  February  1977  event  will  provide  a  unique  opportunity  for 
validation  when  such  events  are  available.  Following  the  tomographic  inversion  and  subsequent  model 
validation,  we  believe  that  the  updated  WINPAK3D  model  will  provide  reliable,  high-accuracy  locations, 
for  small  events  with  sparse  data,  which  will  be  critical  for  CTBT  monitoring  in  the  India-Pakistan  region. 
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